---
title: 'Balancing Disproportionality and Parliament Fragmentation: A Simulation Study of the Mechanical Effects of District Magnitude on Electoral Outcomes'
author: "Micha? Pierzgalski"
date: "June 18, 2018"
output: html_document

#  beamer_presentation: default
#  slidy_presentation: default
---

```{r setup, include=FALSE}

# Create the following directory structure to replicate the analysis:
# |-- sim_replicate_R_Markdown_18_06_18.Rmd
# |-- Figs
# |-- data
#     `-- election_dk_1947.csv
#     `-- clea_1989_dm.csv

## !!! Default random seed used to simulate data = 0

## !!! Depending on the CPU and RAM, computations may take a long time

knitr::opts_chunk$set(fig.width=26, fig.height=12, dpi = 300, fig.path='Figs/', echo=FALSE, warning=F, message=T)

```


# SIMULATION

The following packages are required to conduct the study:

```{r, Load packages, message=FALSE}
library(devtools)

library(ggplot2)
library(lattice)
library(viridis)
library(ggthemes)
library(psych)
library(dplyr)
library(tidyr)
library(truncdist)
library(data.table)
library(cowplot)
library(gtools)

library(stargazer)

# Load disprr - the package that I developed to conduct simulations (the package is still under development (pre-alpha version) and currently it is just the set of interrelated R functions).

# The source code of the package is anonymously hosted at GitHub. Package website at GitHub: https://pierzgal.github.io/disprr/
devtools::install_github("pierzgal/disprr")
library(disprr)

# remove.packages("disprr", lib="~/R/win-library/3.3")
```

----

# Log-normal Distribution Used to Sample Data [logN(mean = 9, sd = 0.8)]

```{r}
options(scipen = 999)

x <- seq(0, 100000, by = 600)
y <-
  dtrunc(
    x,
    spec = "lnorm",
    mean = 9,
    sd = 0.8,
    a = 0,
    b = 100000
  )
xy <- data.frame(Votes = x, Density = y)

xyplot(Density ~ Votes, type = "l", data = xy)


x <- seq(0, 100000, by = 600)
y <-
  dtrunc(
    x,
    spec = "exp",
    rate = 1 / 25000,
    a = 0,
    b = 100000
  )
xy <- data.frame(Votes = x, Density = y)

xyplot(Density ~ Votes, type = "l", data = xy)
```

----

# Empirical Distribution of Vote Shares for Results of "Norway", "Finland", "Switzerland", "Sweden", "Ireland", "Portugal", "Denmark" and "Italy" PR Elections in the period 1991-2016 

```{r}

# WESTERN EUROPE - PR Elections (1991-2016)

clea_1989_dm <-
  read.csv(file = "data/clea_1989_dm.csv", header = TRUE, sep = ",")

ctr_n_1 <- filter(clea_1989_dm, ctr_n %in% c("Portugal", "Denmark", "Italy") )
ctr_n_3 <- filter(clea_1989_dm, ctr_n %in% c("Greece", "Spain") )

ctr_n <- filter(clea_1989_dm, ctr_n %in% c("Norway", "Finland", "Switzerland", "Sweden", "Ireland", "Portugal", "Denmark", "Italy") )

densityplot(~ pvs1, data = ctr_n, xlab = "Vote fractions (constituency level)", ylab = "Density")

densityplot(~ pvs1 | ctr_n, data = ctr_n_1)
densityplot(~ pvs1 | ctr_n, data = ctr_n_2, xlab = "Vote fractions (constituency level)", ylab = "Density")
densityplot(~ pvs1 | ctr_n, data = ctr_n_3)


# DENMARK

DK <-
  read.csv(file = "data/election_dk_1947.csv", header = TRUE, sep = ",")
  
  ggplot( data = subset(DK, yr != 2015 ) ) + geom_density( aes( x = pvs1 ) ) + theme_classic() + xlab("Vote fractions (constituency level) - Denmark : 1947-2011") + ylab("Density")
  
  bwplot( ~ pvs1, breaks = 100, ylab = "Density", xlab = "Vote fractions (constituency level) - Denmark : 1947-2011", jitter = 0.5, alpha = 1, aspect = 0.5, data = subset(DK, yr != 2015 ) )
  
# NP <- DK %>% group_by(yr, cst_n) %>% summarise( np = n() )
# median(NP$np)
  

```

----

# Distributions of Vote Shares and Vote Counts for Simulated Data

```{r}

## For one multi-member district and using the log-normal density function

sample <- sampleElectionData(
  seed = 0,
  dist = "lnorm",
  np = 6,
  nd = 1,
  ne = 500,
  mean = 9,
  sd = 0.8,
#  rate = 1 / 2000,
  max = 100000,
  TS = 20,
  formula_dist = "hh"
  )

vs <- data.frame( Vote_Shares = unlist(sample[[3]]) )
vt <- data.frame( Vote = unlist(sample[['Votes_Total_Party']]) )

histogram( ~ Vote_Shares, breaks = 100, ylab = "Density", xlab = "Vote Shares", aspect = 0.5, data = vs )
bwplot( ~ Vote_Shares, breaks = 100, ylab = "Density", xlab = "Vote Shares", jitter = 0.5, alpha = 1, aspect = 0.5, data = vs )
densityplot( ~ Vote_Shares, breaks = 100, ylab = "Density", xlab = "Vote Shares", aspect = 0.5, data = vs )

histogram( ~ Vote, breaks = 100, ylab = "Density", xlab = "Vote Count", aspect = 0.5, data = vt )
densityplot( ~ Vote, breaks = 100, ylab = "Density", xlab = "Vote Count", aspect = 0.5, data = vt )

```

----

# AGGREGATE-LEVEL MEASURES OF DISPROPORTIONALITY

## GHI ~ district magnitude

```{r, warning=FALSE}
# For 3-party system

sim_aggregate3 <- Disp2(seed = 0, np = 3, dist = "lnorm", mean = 9, sd = 0.8, ne = 200, minTS = 2, maxTS = 20)
```

```{r}
# plot_Disp2(sim_aggregate3)

plot_Disp2(sim_aggregate3, methods = c("DH", "SL", "H", "Danish", "MSL", "Imperiali"))[["plot_GHI"]]
ggsave("sim_aggregate3_1.pdf", width = 30, height = 20, units = "cm")

plot_Disp2(sim_aggregate3, methods = c("DH", "SL", "H", "Danish", "MSL", "Imperiali"))[["plot_NPP"]]
ggsave("sim_aggregate3_npp.pdf", width = 30, height = 19, units = "cm")

plot_Disp2(sim_aggregate3, methods = c("DH", "SL", "H", "Danish", "MSL", "Imperiali"))[["plot_ENPP"]]
ggsave("sim_aggregate3_enpp.pdf", width = 34, height = 19, units = "cm")
```


```{r, warning=FALSE}
# For 6-party system

sim_aggregate6 <- Disp2(np = 6, dist = "lnorm", mean = 9, sd = 0.8, ne = 200, minTS = 2, maxTS = 20)
```

```{r}
# plot_Disp2(sim_aggregate6)

plot_Disp2(sim_aggregate6, methods = c("DH", "SL", "H", "Danish", "MSL", "Imperiali"))[["plot_GHI"]]
ggsave("sim_aggregate6_1.pdf", width = 30, height = 20, units = "cm")

plot_Disp2(sim_aggregate6, methods = c("DH", "SL", "H", "Danish", "MSL", "Imperiali"))[["plot_NPP"]]
ggsave("sim_aggregate6_npp.pdf", width = 30, height = 19, units = "cm")

plot_Disp2(sim_aggregate6, methods = c("DH", "SL", "H", "Danish", "MSL", "Imperiali"))[["plot_ENPP"]]
ggsave("sim_aggregate6_enpp.pdf", width = 34, height = 19, units = "cm")
```


```{r, warning=FALSE}
# For 9-party system

sim_aggregate9 <- Disp2(np = 9, dist = "lnorm", mean = 9, sd = 0.8, ne = 200, minTS = 2, maxTS = 20)
```

```{r}
# plot_Disp2(sim_aggregate9)

plot_Disp2(sim_aggregate9, methods = c("DH", "SL", "H", "Danish", "MSL", "Imperiali"))[["plot_GHI"]]
ggsave("sim_aggregate9_1.pdf", width = 30, height = 20, units = "cm")

plot_Disp2(sim_aggregate9, methods = c("DH", "SL", "H", "Danish", "MSL", "Imperiali"))[["plot_NPP"]]
ggsave("sim_aggregate9_npp.pdf", width = 30, height = 19, units = "cm")

plot_Disp2(sim_aggregate9, methods = c("DH", "SL", "H", "Danish", "MSL", "Imperiali"))[["plot_ENPP"]]
ggsave("sim_aggregate9_enpp.pdf", width = 34, height = 19, units = "cm")
```


```{r, warning=FALSE}
# For 12-party system

sim_aggregate12 <- Disp2(np = 12, dist = "lnorm", mean = 9, sd = 0.8, ne = 100, minTS = 2, maxTS = 20)
```

```{r}
plot_Disp2(sim_aggregate12)

plot_Disp2(sim_aggregate12, methods = c("DH", "SL", "H", "Danish", "MSL", "Imperiali"))[["plot_GHI"]]
ggsave("sim_aggregate12_1.pdf", width = 30, height = 20, units = "cm")

plot_Disp2(sim_aggregate12, methods = c("DH", "SL", "H", "Danish", "MSL", "Imperiali"))[["plot_NPP"]]
ggsave("sim_aggregate12_npp.pdf", width = 30, height = 19, units = "cm")

plot_Disp2(sim_aggregate12, methods = c("DH", "SL", "H", "Danish", "MSL", "Imperiali"))[["plot_ENPP"]]
ggsave("sim_aggregate12_enpp.pdf", width = 34, height = 19, units = "cm")
```

----

### 3-party system

```{r}
sim_aggregate3[["summary"]]
```

----

### 6-party system

```{r}
sim_aggregate6[["summary"]]

```

----

### 9-party system

```{r}
sim_aggregate9[["summary"]]

```

----

### 12-party system

```{r}
sim_aggregate12[["summary"]]

```

----

# Model (log-level with interactions using pooled data)

```{r}

sim_aggregate_all <- dplyr::bind_rows(sim_aggregate3[["summary"]], sim_aggregate6[["summary"]], sim_aggregate9[["summary"]], sim_aggregate12[["summary"]])

sim_aggregate_all <- dplyr::bind_rows(sim_aggregate9[["summary"]], sim_aggregate3[["summary"]]) # ...


sim_aggregate_all <- within(sim_aggregate_all, method <- relevel(factor(method), ref = "H"))


model_all <- lm( log(GHI) ~ DM*factor(method) + DM*factor(NP), data = filter(sim_aggregate_all) )
summary(model_all)
  
# The results

# Call:
# lm(formula = log(GHI) ~ DM * factor(method) + DM * factor(NP), 
#     data = filter(sim_aggregate_all))
# 
# Residuals:
#     Min      1Q  Median      3Q     Max 
# -4.0954 -0.1779  0.0092  0.2146  1.3742 
# 
# Coefficients:
#                              Estimate Std. Error  t value Pr(>|t|)    
# (Intercept)                -2.0359772  0.0051425 -395.913  < 2e-16 ***
# DM                         -0.1221038  0.0004185 -291.772  < 2e-16 ***
# factor(method)Danish        0.0171552  0.0059380    2.889 0.003865 ** 
# factor(method)DH            0.1752764  0.0059380   29.518  < 2e-16 ***
# factor(method)Imperiali     0.4133517  0.0059380   69.611  < 2e-16 ***
# factor(method)MSL           0.1375372  0.0059380   23.162  < 2e-16 ***
# factor(method)SL            0.0207441  0.0059380    3.493 0.000477 ***
# factor(NP)6                 0.4757755  0.0048484   98.131  < 2e-16 ***
# factor(NP)9                 0.6460878  0.0048484  133.258  < 2e-16 ***
# factor(NP)12                0.7412916  0.0048484  152.894  < 2e-16 ***
# DM:factor(method)Danish     0.0040976  0.0004832    8.480  < 2e-16 ***
# DM:factor(method)DH         0.0117936  0.0004832   24.406  < 2e-16 ***
# DM:factor(method)Imperiali  0.0407341  0.0004832   84.295  < 2e-16 ***
# DM:factor(method)MSL       -0.0027553  0.0004832   -5.702 1.19e-08 ***
# DM:factor(method)SL        -0.0002258  0.0004832   -0.467 0.640376    
# DM:factor(NP)6              0.0037491  0.0003946    9.502  < 2e-16 ***
# DM:factor(NP)9              0.0089002  0.0003946   22.557  < 2e-16 ***
# DM:factor(NP)12             0.0132816  0.0003946   33.662  < 2e-16 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.3648 on 227982 degrees of freedom
# Multiple R-squared:  0.8061,	Adjusted R-squared:  0.8061 
# F-statistic: 5.576e+04 on 17 and 227982 DF,  p-value: < 2.2e-16

# ----

# LaTeX table
stargazer(model_all, style = "apsr", single.row = T, ci = T)


### (TEST)

# # method : H, NP = 3, DM = 8
# 
# GHI = exp( Intercept + cDM * 8 + cDH * 0 + cNP6 * 0 + cDM_DH * 8 * 0 + cDM_NP6 * 8 * 0 )
# GHI
# 
# #
# 
# # method : H -> DH | NP = 6, DM = 6
# 
# GHI0 = exp( Intercept + cDM * 6 + cDH * 0 + cNP6 * 1 + cDM_DH * 6 * 0 + cDM_NP6 * 6 * 1 )
# GHI1 = exp( Intercept + cDM * 6 + cDH * 1 + cNP6 * 1 + cDM_DH * 6 * 1 + cDM_NP6 * 6 * 1 )
# 
# deltaGHI = (GHI1-GHI0) / GHI0
# deltaGHI
# 
# exp(cDH + cDM_DH * 6 * 1 ) - 1
# 
# 
# ###
# 
# # DM : 3 -> 6 | NP = 3, method = H
# 
# GHI0 = exp( Intercept + cDM * 3 + cDH * 0 + cNP6 * 0 + cDM_DH * 3 * 0 + cDM_NP6 * 3 * 0 )
# GHI1 = exp( Intercept + cDM * 6 + cDH * 0 + cNP6 * 0 + cDM_DH * 6 * 0 + cDM_NP6 * 6 * 0 )
# 
# deltaGHI = (GHI1-GHI0) / GHI0
# deltaGHI
# 
# exp(cDM * 3 + cDM_DH * 3 * 0 + cDM_NP6 * 3 * 0) - 1
# 
# # DM : 6 -> 9 | NP = 6, method = DH
# 
# GHI0 = exp( Intercept + cDM * 6 + cDH * 1 + cNP6 * 1 + cDM_DH * 6 * 1 + cDM_NP6 * 6 * 1 )
# GHI1 = exp( Intercept + cDM * 9 + cDH * 1 + cNP6 * 1 + cDM_DH * 9 * 1 + cDM_NP6 * 9 * 1 )
# 
# deltaGHI = (GHI1-GHI0) / GHI0
# deltaGHI
# 
# exp(cDM * 3 + cDM_DH * 3 * 1 + cDM_NP6 * 3 * 1) - 1
# 
# # DM : 9 -> 6
# 
# GHI0 = exp( Intercept + cDM * 6 + cDH * 1 + cNP6 * 1 + cDM_DH * 6 * 1 + cDM_NP6 * 6 * 1 )
# GHI1 = exp( Intercept + cDM * 9 + cDH * 1 + cNP6 * 1 + cDM_DH * 9 * 1 + cDM_NP6 * 9 * 1 )
# 
# deltaGHI = (GHI0-GHI1) / GHI1
# deltaGHI
# 
# exp( -(cDM * 3 + cDM_DH * 3 * 1 + cDM_NP6 * 3 * 1) ) - 1

### END


# Plot log-level and nls
  
## Compare models (GHI ~ DM | without moderators) 
  
model_all2 <- lm( log(GHI) ~ DM, data = filter(sim_aggregate_all) )
summary(model_all2)

model_all3 <-
    nls(
      GHI ~ C * exp(alpha * DM),
      start = list(C = 0.3, alpha = -0.1),
      data = filter(sim_aggregate_all)
    )

all2 = filter(sim_aggregate_all)
all2 = dplyr::mutate(all2, GHI_predicted = predict(model_all3), GHI_loglevel = exp( predict( model_all2 ) ) ) # Compute predicted values based on the models

  plot_GHI_all <-
    ggplot2::ggplot(data = all2) + ggplot2::geom_boxplot(
      ggplot2::aes(x = as.factor(DM),
                   y = GHI),
      lwd = 0.25,
      fatten = 0.4,
      outlier.size = 0.6
    ) + ggplot2::geom_line(
      ggplot2::aes(x = as.factor(DM), y = GHI_predicted, group = 1),
      size = 0.35,
      colour = "blue"
    ) + ggplot2::geom_line(
      ggplot2::aes(x = as.factor(DM), y = GHI_loglevel, group = 1),
      size = 0.35,
      colour = "red",
      linetype = "dashed"
    ) + xlab("DM")

  plot_GHI_all

```

----

## Interpreting the log-level model

```{r}

# Plot relative effects of DM changes (deltaDM = 1) for GHI : for various configurations of apportionment methods

model_ll_h <- lm( log(GHI) ~ DM, data = filter(sim_aggregate_all, method == "H") )
model_ll_dh <- lm( log(GHI) ~ DM, data = filter(sim_aggregate_all, method == "DH") )
model_ll_imp <- lm( log(GHI) ~ DM, data = filter(sim_aggregate_all, method == "Imperiali") )
model_ll_sl <- lm( log(GHI) ~ DM, data = filter(sim_aggregate_all, method == "SL") )

ll_h = filter(sim_aggregate_all, method == "H")
ll_dh = filter(sim_aggregate_all, method == "DH")
ll_imp = filter(sim_aggregate_all, method == "Imperiali")
ll_sl = filter(sim_aggregate_all, method == "SL")

ll_h = dplyr::mutate(ll_h, GHI_loglevel = exp( predict( model_ll_h ) ) )
ll_dh = dplyr::mutate(ll_dh, GHI_loglevel = exp( predict( model_ll_dh ) ) )
ll_imp = dplyr::mutate(ll_imp, GHI_loglevel = exp( predict( model_ll_imp ) ) )
ll_sl = dplyr::mutate(ll_sl, GHI_loglevel = exp( predict( model_ll_sl ) ) )

ll_all <- bind_rows(ll_h, ll_dh, ll_imp, ll_sl)

ll_all <- filter( ll_all, DM >= 4, method != "H" ) # ...

ggplot(data = ll_all) + geom_boxplot(
  ggplot2::aes(x = as.factor(DM),
  y = GHI, fill = method),
  lwd = 0.25,
  fatten = 0.4,
  outlier.size = 0.6
  ) + xlab("DM") + geom_hline(
  yintercept = c(0, 0.05, 0.1),
  size = 0.35,
  linetype = "longdash",
  colour = "blue") + facet_wrap(~ NP)

ggplot(data = ll_all) + geom_text(aes(
  x = factor(DM),
  y = GHI_loglevel,
  label = signif(GHI_loglevel, digits = 3),
  color = method
  ), check_overlap = F) + xlab("DM") + geom_hline(
  yintercept = c(0, 0.05, 0.1, 0.15, 0.2),
  size = 0.35,
  linetype = "longdash",
  colour = "blue"
  ) + ylab("Gallagher index") + scale_fill_viridis(discrete = T)
  
```

---

### How an allocation formula and a party system format moderate the impact of DM on disproportionality

```{r}

## Define vars containg regression coefficient values

Intercept =                          -2.0359772
cDM =                                   -0.1221038
cDanish =                  0.0171552
cDH =                      0.1752764
cImperiali =               0.4133517
cMSL =                     0.1375372
cSL =                      0.0207441
cNP6 =                           0.4757755
cNP9 =                           0.6460878
cNP12 =                          0.7412916
cDM_Danish =               0.0040976
cDM_DH =                   0.0117936
cDM_Imperiali =            0.0407341
cDM_MSL =                 -0.0027553
cDM_SL =                  -0.0002258
cDM_NP6 =                        0.0037491
cDM_NP9 =                        0.0089002
cDM_NP12 =                       0.0132816

```

----

```{r}

### DM (-) | method

#### 1

jump = 3 # define DM change (decrease DM by ... seats)


##### method = H

DM_contrast3h <- data.frame( method = "H", NP = 3, decrease = exp(-(cDM * jump + cDM_Danish * jump * 0 + cDM_NP6 * jump * 0)) - 1 ) # NP = 3
DM_contrast6h <- data.frame( method = "H", NP = 6, decrease = exp(-(cDM * jump + cDM_Danish * jump * 0 + cDM_NP6 * jump * 1)) - 1 )
DM_contrast9h <- data.frame( method = "H", NP = 9, decrease = exp(-(cDM * jump + cDM_Danish * jump * 0 + cDM_NP9 * jump * 1)) - 1 )
DM_contrast12h <- data.frame( method = "H", NP = 12, decrease = exp(-(cDM * jump + cDM_Danish * jump * 0 + cDM_NP12 * jump * 1 )) - 1 )

##### method = Danish

DM_contrast3d <- data.frame( method = "Danish", NP = 3, decrease = exp(-(cDM * jump + cDM_Danish * jump * 1 + cDM_NP6 * jump * 0)) - 1 ) # NP = 3
DM_contrast6d <- data.frame( method = "Danish", NP = 6, decrease = exp(-(cDM * jump + cDM_Danish * jump * 1 + cDM_NP6 * jump * 1)) - 1 )
DM_contrast9d <- data.frame( method = "Danish", NP = 9, decrease = exp(-(cDM * jump + cDM_Danish * jump * 1 + cDM_NP9 * jump * 1)) - 1 )
DM_contrast12d <- data.frame( method = "Danish", NP = 12, decrease = exp(-(cDM * jump + cDM_Danish * jump * 1 + cDM_NP12 * jump * 1 )) - 1 )

##### method = DH

DM_contrast3dh <- data.frame( method = "DH", NP = 3, decrease = exp(-(cDM * jump + cDM_DH * jump * 1 + cDM_NP6 * jump * 0 )) - 1 ) # NP = 3
DM_contrast6dh <- data.frame( method = "DH", NP = 6, decrease = exp(-(cDM * jump + cDM_DH * jump * 1 + cDM_NP6 * jump * 1 )) - 1 ) 
DM_contrast9dh <- data.frame( method = "DH", NP = 9, decrease = exp(-(cDM * jump + cDM_DH * jump * 1 + cDM_NP9 * jump * 1 )) - 1 )
DM_contrast12dh <- data.frame( method = "DH", NP = 12, decrease = exp(-(cDM * jump + cDM_DH * jump * 1 + cDM_NP12 * jump * 1 )) - 1 )

##### method = SL

DM_contrast3sl <- data.frame( method = "SL", NP = 3, decrease = exp(-(cDM * jump + cDM_SL * jump * 1 + cDM_NP6 * jump * 0)) - 1 ) # NP = 3
DM_contrast6sl <- data.frame( method = "SL", NP = 6, decrease = exp(-(cDM * jump + cDM_SL * jump * 1 + cDM_NP6 * jump * 1)) - 1 )
DM_contrast9sl <- data.frame( method = "SL", NP = 9, decrease = exp(-(cDM * jump + cDM_SL * jump * 1 + cDM_NP9 * jump * 1)) - 1 )
DM_contrast12sl <- data.frame( method = "SL", NP = 12, decrease = exp(-(cDM * jump + cDM_SL * jump * 1 + cDM_NP12 * jump * 1 )) - 1 )

##### method = MSL

DM_contrast3msl <- data.frame( method = "MSL", NP = 3, decrease = exp(-(cDM * jump + cDM_MSL * jump * 1 + cDM_NP6 * jump * 0)) - 1 ) # NP = 3
DM_contrast6msl <- data.frame( method = "MSL", NP = 6, decrease = exp(-(cDM * jump + cDM_MSL * jump * 1 + cDM_NP6 * jump * 1)) - 1 )
DM_contrast9msl <- data.frame( method = "MSL", NP = 9, decrease = exp(-(cDM * jump + cDM_MSL * jump * 1 + cDM_NP9 * jump * 1)) - 1 )
DM_contrast12msl <- data.frame( method = "MSL", NP = 12, decrease = exp(-(cDM * jump + cDM_MSL * jump * 1 + cDM_NP12 * jump * 1 )) - 1 )

##### method = Imperiali

DM_contrast3i <- data.frame( method = "Imperiali", NP = 3, decrease = exp(-(cDM * jump + cDM_Imperiali * jump * 1 + cDM_NP6 * jump * 0)) - 1 ) # NP = 3
DM_contrast6i <- data.frame( method = "Imperiali", NP = 6, decrease = exp(-(cDM * jump + cDM_Imperiali * jump * 1 + cDM_NP6 * jump * 1)) - 1 )
DM_contrast9i <- data.frame( method = "Imperiali", NP = 9, decrease = exp(-(cDM * jump + cDM_Imperiali * jump * 1 + cDM_NP9 * jump * 1)) - 1 )
DM_contrast12i <- data.frame( method = "Imperiali", NP = 12, decrease = exp(-(cDM * jump + cDM_Imperiali * jump * 1 + cDM_NP12 * jump * 1 )) - 1 )


# ---

DM_contrast <- bind_rows(DM_contrast3d, DM_contrast6d, DM_contrast9d, DM_contrast12d, DM_contrast3dh, DM_contrast6dh, DM_contrast9dh, DM_contrast12dh, DM_contrast3sl, DM_contrast6sl, DM_contrast9sl, DM_contrast12sl, DM_contrast3h, DM_contrast6h, DM_contrast9h, DM_contrast12h, DM_contrast3msl, DM_contrast6msl, DM_contrast9msl, DM_contrast12msl, DM_contrast3i, DM_contrast6i, DM_contrast9i, DM_contrast12i)

# Plot
ggplot(data = filter(DM_contrast), aes(y = factor(NP), x = method)) + geom_tile(aes(fill = decrease)) + geom_text(aes(label = signif(decrease, 4))) + scale_fill_gradient("Relative increase\nin GHI | DeltaDM = -3", low = "white", high = "grey") + ylab("NP") + xlab("Method")

```

----

```{r}

### DM (+) | method

#### 1

jump = 3 # define DM change (increse DM by ... seats)


##### method = H

DM_contrast3h <- data.frame( method = "H", NP = 3, decrease = exp(cDM * jump + cDM_Danish * jump * 0 + cDM_NP6 * jump * 0) - 1 ) # NP = 3
DM_contrast6h <- data.frame( method = "H", NP = 6, decrease = exp(cDM * jump + cDM_Danish * jump * 0 + cDM_NP6 * jump * 1) - 1 )
DM_contrast9h <- data.frame( method = "H", NP = 9, decrease = exp(cDM * jump + cDM_Danish * jump * 0 + cDM_NP9 * jump * 1) - 1 )
DM_contrast12h <- data.frame( method = "H", NP = 12, decrease = exp(cDM * jump + cDM_Danish * jump * 0 + cDM_NP12 * jump * 1 ) - 1 )

##### method = Danish

DM_contrast3d <- data.frame( method = "Danish", NP = 3, decrease = exp(cDM * jump + cDM_Danish * jump * 1 + cDM_NP6 * jump * 0) - 1 ) # NP = 3
DM_contrast6d <- data.frame( method = "Danish", NP = 6, decrease = exp(cDM * jump + cDM_Danish * jump * 1 + cDM_NP6 * jump * 1) - 1 )
DM_contrast9d <- data.frame( method = "Danish", NP = 9, decrease = exp(cDM * jump + cDM_Danish * jump * 1 + cDM_NP9 * jump * 1) - 1 )
DM_contrast12d <- data.frame( method = "Danish", NP = 12, decrease = exp(cDM * jump + cDM_Danish * jump * 1 + cDM_NP12 * jump * 1 ) - 1 )

##### method = DH

DM_contrast3dh <- data.frame( method = "DH", NP = 3, decrease = exp(cDM * jump + cDM_DH * jump * 1 + cDM_NP6 * jump * 0 ) - 1 ) # NP = 3
DM_contrast6dh <- data.frame( method = "DH", NP = 6, decrease = exp(cDM * jump + cDM_DH * jump * 1 + cDM_NP6 * jump * 1 ) - 1 )
DM_contrast9dh <- data.frame( method = "DH", NP = 9, decrease = exp(cDM * jump + cDM_DH * jump * 1 + cDM_NP9 * jump * 1 ) - 1 )
DM_contrast12dh <- data.frame( method = "DH", NP = 12, decrease = exp(cDM * jump + cDM_DH * jump * 1 + cDM_NP12 * jump * 1 ) - 1 )

##### method = SL

DM_contrast3sl <- data.frame( method = "SL", NP = 3, decrease = exp(cDM * jump + cDM_SL * jump * 1 + cDM_NP6 * jump * 0) - 1 ) # NP = 3
DM_contrast6sl <- data.frame( method = "SL", NP = 6, decrease = exp(cDM * jump + cDM_SL * jump * 1 + cDM_NP6 * jump * 1) - 1 )
DM_contrast9sl <- data.frame( method = "SL", NP = 9, decrease = exp(cDM * jump + cDM_SL * jump * 1 + cDM_NP9 * jump * 1) - 1 )
DM_contrast12sl <- data.frame( method = "SL", NP = 12, decrease = exp(cDM * jump + cDM_SL * jump * 1 + cDM_NP12 * jump * 1 ) - 1 )

##### method = MSL

DM_contrast3msl <- data.frame( method = "MSL", NP = 3, decrease = exp(cDM * jump + cDM_MSL * jump * 1 + cDM_NP6 * jump * 0) - 1 ) # NP = 3
DM_contrast6msl <- data.frame( method = "MSL", NP = 6, decrease = exp(cDM * jump + cDM_MSL * jump * 1 + cDM_NP6 * jump * 1) - 1 )
DM_contrast9msl <- data.frame( method = "MSL", NP = 9, decrease = exp(cDM * jump + cDM_MSL * jump * 1 + cDM_NP9 * jump * 1) - 1 )
DM_contrast12msl <- data.frame( method = "MSL", NP = 12, decrease = exp(cDM * jump + cDM_MSL * jump * 1 + cDM_NP12 * jump * 1 ) - 1 )

##### method = Imperiali

DM_contrast3i <- data.frame( method = "Imperiali", NP = 3, decrease = exp(cDM * jump + cDM_Imperiali * jump * 1 + cDM_NP6 * jump * 0) - 1 ) # NP = 3
DM_contrast6i <- data.frame( method = "Imperiali", NP = 6, decrease = exp(cDM * jump + cDM_Imperiali * jump * 1 + cDM_NP6 * jump * 1) - 1 )
DM_contrast9i <- data.frame( method = "Imperiali", NP = 9, decrease = exp(cDM * jump + cDM_Imperiali * jump * 1 + cDM_NP9 * jump * 1) - 1 )
DM_contrast12i <- data.frame( method = "Imperiali", NP = 12, decrease = exp(cDM * jump + cDM_Imperiali * jump * 1 + cDM_NP12 * jump * 1 ) - 1 )


# ---

DM_contrast <-
  bind_rows(
  DM_contrast3d,
  DM_contrast6d,
  DM_contrast9d,
  DM_contrast12d,
  DM_contrast3dh,
  DM_contrast6dh,
  DM_contrast9dh,
  DM_contrast12dh,
  DM_contrast3sl,
  DM_contrast6sl,
  DM_contrast9sl,
  DM_contrast12sl,
  DM_contrast3h,
  DM_contrast6h,
  DM_contrast9h,
  DM_contrast12h,
  DM_contrast3msl,
  DM_contrast6msl,
  DM_contrast9msl,
  DM_contrast12msl,
  DM_contrast3i,
  DM_contrast6i,
  DM_contrast9i,
  DM_contrast12i
  )

# Plot
ggplot(data = filter(DM_contrast), aes(y = factor(NP), x = method)) + geom_tile(aes(fill = decrease)) + geom_text(aes(label = signif(decrease, 4))) + scale_fill_gradient("Relative decrease\nin GHI | DeltaDM = +3", low = "white", high = "grey") + ylab("NP") + xlab("Method")

```

----

```{r}

### H -> DH, ... | DM

DH_contrast <- data.frame( DM = c(2:20), DH = exp(cDH + cDM_DH * c(2:20) * 1 ) - 1 ) # ADD cDM_NP6 * 0
Danish_contrast <- data.frame( Danish = exp(cDanish + cDM_Danish * c(2:20) * 1 ) - 1 )
Imperiali_contrast <- data.frame( Imperiali = exp(cImperiali + cDM_Imperiali * c(2:20) * 1 ) - 1 )
MSL_contrast <- data.frame( MSL = exp(cMSL + cDM_MSL * c(2:20) * 1 ) - 1 )
SL_contrast <- data.frame( SL = exp(cSL + cDM_SL * c(2:20) * 1 ) - 1 )

contrasts = bind_cols(DH_contrast, Danish_contrast, Imperiali_contrast, MSL_contrast, SL_contrast)
contrasts_long <- gather(contrasts, formula, contrast, DH:SL, factor_key=TRUE)

mean(DH_contrast$DH)


### Plot

ggplot(data = filter(contrasts_long), aes(x = factor(DM), y = formula)) + geom_tile(aes(fill = contrast)) + geom_text(aes(label = round(contrast, 2))) + scale_fill_gradient("Relative change | Hamilton", low = "white", high = "grey") + theme_classic() + xlab("DM") + ylab("Method")

```

----

```{r}

### NP = 3 -> 6, 9, 12 | DM

contrast6 <- data.frame( DM = c(2:20), NP6 = exp(cNP6 + cDM_NP6 * c(2:20) * 1 ) - 1 ) # ADD cDM_NP6 * 0
contrast9 <- data.frame( NP9 = exp(cNP9 + cDM_NP9 * c(2:20) * 1 ) - 1 )
contrast12 <- data.frame( NP12 = exp(cNP12 + cDM_NP12 * c(2:20) * 1 ) - 1 )

contrasts = bind_cols(contrast6, contrast9, contrast12)
contrasts_long <- gather(contrasts, formula, contrast, NP6:NP12, factor_key=TRUE)


### Plot

ggplot(data = filter(contrasts_long), aes(x = factor(DM), y = formula)) + geom_tile(aes(fill = contrast)) + geom_text(aes(label = round(contrast, 2))) + scale_fill_gradient("Relative change | NP = 3", low = "white", high = "grey") + theme_classic() + xlab("DM") + ylab("# of Parties")

```

----

# Models (non-linear : Disp ~ DM | ...)

## Model - d'Hondt

```{r}
# dHondt

model_dh3 <- summary(sim_aggregate3[["Model_DH"]])
model_dh3

# x = seq(3, 20, by = 1)
# y_dh = model_dh3$coefficients[1] * exp(model_dh3$coefficients[2]*x)
# 
# dydx_dh = model_dh3$coefficients[1] * model_dh3$coefficients[2] * exp(model_dh3$coefficients[2]*x)
# Elasticity = model_dh3$coefficients[2] * x
# 
# scatterplot_dh3 <- data.frame(GHI_predicted = y_dh, DM = x, dy_dx = dydx_dh, Elasticity, method = "DH")
# lattice::xyplot(Elasticity ~ DM, data = scatterplot_dh3)
# lattice::xyplot(GHI_predicted ~ factor(DM), data = scatterplot_dh3)

```

----

```{r}
# dHondt

model_dh6 <- summary(sim_aggregate6[["Model_DH"]])
model_dh6

# x = seq(3, 20, by = 1)
# y_dh = model_dh6$coefficients[1] * exp(model_dh6$coefficients[2]*x)
# 
# dydx_dh = model_dh6$coefficients[1] * model_dh6$coefficients[2] * exp(model_dh6$coefficients[2]*x)
# 
# scatterplot_dh6 <- data.frame(GHI_predicted = y_dh, DM = x, dy_dx = dydx_dh, method = "DH")
# lattice::xyplot(dy_dx ~ DM, data = scatterplot_dh6)
# lattice::xyplot(GHI_predicted ~ factor(DM), xlab = 'DM', data = scatterplot_dh6)

```

----

```{r}
# dHondt

model_dh9 <- summary(sim_aggregate9[["Model_DH"]])
model_dh9

# x = seq(3, 20, by = 1)
# y_dh9 = model_dh9$coefficients[1] * exp(model_dh9$coefficients[2]*x)
# 
# dydx_dh9 = model_dh9$coefficients[1] * model_dh9$coefficients[2] * exp(model_dh9$coefficients[2]*x)
# 
# scatterplot_dh9 <- data.frame(GHI_predicted = y_dh9, DM = x, dy_dx = dydx_dh9, method = "DH")
# lattice::xyplot(dy_dx ~ DM, data = scatterplot_dh9)
# lattice::xyplot(GHI_predicted ~ factor(DM), data = scatterplot_dh9)

```

----

```{r}
# dHondt

model_dh12 <- summary(sim_aggregate12[["Model_DH"]])
model_dh12

# x = seq(3, 20, by = 1)
# y_dh12 = model_dh12$coefficients[1] * exp(model_dh12$coefficients[2]*x)
# 
# dydx_dh12 = model_dh12$coefficients[1] * model_dh12$coefficients[2] * exp(model_dh12$coefficients[2]*x)
# 
# scatterplot_dh12 <- data.frame(GHI_predicted = y_dh12, DM = x, dy_dx = dydx_dh12, method = "DH", np = 12)
# lattice::xyplot(dy_dx ~ DM, data = scatterplot_dh12)
# lattice::xyplot(GHI_predicted ~ factor(DM), data = scatterplot_dh12)

```

----

## Model - Sainte-Lague

```{r}
# SL

model_sl3 <- summary(sim_aggregate3[["Model_SL"]])
model_sl3

# x = seq(3, 20, by = 1)
# y_sl = exp(model_sl3$coefficients[1] + model_sl3$coefficients[2]*x)
# 
# dydx_sl = model_sl3$coefficients[1] * exp(model_sl3$coefficients[1] + model_sl3$coefficients[2]*x)
# 
# scatterplot_sl3 <- data.frame(GHI_predicted = y_sl, DM = x, dy_dx = dydx_sl, method = "SL")
# lattice::xyplot(dy_dx ~ DM, data = scatterplot_sl3)
# lattice::xyplot(GHI_predicted ~ factor(DM), data = scatterplot_sl3)

```

----

```{r}
# SL

model_sl6 <- summary(sim_aggregate6[["Model_SL"]])
model_sl6

# x = seq(3, 20, by = 1)
# y_sl = exp(model_sl6$coefficients[1] + model_sl6$coefficients[2]*x)
# 
# dydx_sl = model_sl6$coefficients[1] * exp(model_sl6$coefficients[1] + model_sl6$coefficients[2]*x)
# 
# scatterplot_sl6 <- data.frame(GHI_predicted = y_sl, DM = x, dy_dx = dydx_sl, method = "SL")
# lattice::xyplot(dy_dx ~ DM, data = scatterplot_sl6)
# lattice::xyplot(GHI_predicted ~ factor(DM), data = scatterplot_sl6)

```

----

```{r}
# SL

model_sl9 <- summary(sim_aggregate9[["Model_SL"]])
model_sl9

# x = seq(3, 20, by = 1)
# y_sl = exp(model_sl9$coefficients[1] + model_sl9$coefficients[2]*x)
# 
# dydx_sl = model_sl9$coefficients[1] * exp(model_sl9$coefficients[1] + model_sl9$coefficients[2]*x)
# 
# scatterplot_sl9 <- data.frame(GHI_predicted = y_sl, DM = x, dy_dx = dydx_sl, method = "SL")
# lattice::xyplot(dy_dx ~ DM, data = scatterplot_sl9)
# lattice::xyplot(GHI_predicted ~ factor(DM), data = scatterplot_sl9)

```

----

```{r}
# SL

model_sl12 <- summary(sim_aggregate12[["Model_SL"]])
model_sl12

# x = seq(3, 20, by = 1)
# y_sl12 = exp(model_sl12$coefficients[1] + model_sl12$coefficients[2]*x)
# 
# dydx_sl12 = model_sl12$coefficients[1] * exp(model_sl12$coefficients[1] + model_sl12$coefficients[2]*x)
# 
# scatterplot_sl12 <- data.frame(GHI_predicted = y_sl12, DM = x, dy_dx = dydx_sl12, method = "SL", np = 12)
# lattice::xyplot(dy_dx ~ DM, data = scatterplot_sl12)
# lattice::xyplot(GHI_predicted ~ factor(DM), data = scatterplot_sl12)

```

----

## Model - Imperiali

```{r}
# Imperiali

model_imp3 <- summary(sim_aggregate3[["Model_Imperiali"]])
model_imp3

# x = seq(3, 20, by = 1)
# 
# y_imp3 = model_imp3$coefficients[1] * exp(model_imp3$coefficients[2]*x)
# 
# dydx_imp3 = model_imp3$coefficients[1] * model_imp3$coefficients[2] * exp(model_imp3$coefficients[2]*x)
# 
# scatterplot_imp3 <- data.frame(GHI_predicted = y_imp3, DM = x, dy_dx = dydx_imp3, method = "Imperiali", np = 3)
# lattice::xyplot(dy_dx ~ DM, data = scatterplot_imp3)
# lattice::xyplot(GHI_predicted ~ factor(DM), data = scatterplot_imp3)

```

----

```{r}
# Imperiali

model_imp6 <- summary(sim_aggregate6[["Model_Imperiali"]])
model_imp6

# x = seq(3, 20, by = 1)
# y_imp6 = model_imp6$coefficients[1] * exp(model_imp6$coefficients[2]*x)
# 
# dydx_imp6 = model_imp6$coefficients[1] * model_imp6$coefficients[2] * exp(model_imp6$coefficients[2]*x)
# 
# scatterplot_imp6 <- data.frame(GHI_predicted = y_imp6, DM = x, dy_dx = dydx_imp6, method = "Imperiali", np = 6)
# lattice::xyplot(dy_dx ~ DM, data = scatterplot_imp6)
# lattice::xyplot(GHI_predicted ~ factor(DM), data = scatterplot_imp6)

```

----

```{r}
# Imperiali

model_imp9 <- summary(sim_aggregate9[["Model_Imperiali"]])
model_imp9

# x = seq(3, 20, by = 1)
# #y_imp = exp(model_imp$coefficients[1] + model_imp$coefficients[2]*x)
# y_imp9 = model_imp9$coefficients[1] * exp(model_imp9$coefficients[2]*x)
# 
# dydx_imp9 = model_imp9$coefficients[1] * model_imp9$coefficients[2] * exp(model_imp9$coefficients[2]*x)
# 
# scatterplot_imp9 <- data.frame(GHI_predicted = y_imp9, DM = x, dy_dx = dydx_imp9, method = "Imperiali", np = 9)
# lattice::xyplot(dy_dx ~ DM, data = scatterplot_imp9)
# lattice::xyplot(GHI_predicted ~ factor(DM), data = scatterplot_imp9)

```

----

```{r}
# Imperiali

model_imp12 <- summary(sim_aggregate12[["Model_Imperiali"]])
model_imp12

# x = seq(3, 20, by = 1)
# 
# y_imp12 = model_imp12$coefficients[1] * exp(model_imp12$coefficients[2]*x)
# 
# dydx_imp12 = model_imp12$coefficients[1] * model_imp12$coefficients[2] * exp(model_imp12$coefficients[2]*x)
# 
# scatterplot_imp12 <- data.frame(GHI_predicted = y_imp12, DM = x, dy_dx = dydx_imp12, method = "Imperiali", np = 12)
# lattice::xyplot(dy_dx ~ DM, data = scatterplot_imp12)
# lattice::xyplot(GHI_predicted ~ factor(DM), data = scatterplot_imp12)

```

---

```{r}

# Danish

model_Danish3 <- summary(sim_aggregate3[["Model_Danish"]])
model_Danish3

model_Danish6 <- summary(sim_aggregate6[["Model_Danish"]])
model_Danish6

model_Danish9 <- summary(sim_aggregate9[["Model_Danish"]])
model_Danish9

model_Danish12 <- summary(sim_aggregate12[["Model_Danish"]])
model_Danish12

# H

model_h3 <- summary(sim_aggregate3[["Model_Hamilton"]])
model_h3

model_h6 <- summary(sim_aggregate6[["Model_Hamilton"]])
model_h6

model_h9 <- summary(sim_aggregate9[["Model_Hamilton"]])
model_h9

model_h12 <- summary(sim_aggregate12[["Model_Hamilton"]])
model_h12

# MSL

model_msl3 <- summary(sim_aggregate3[["Model_MSL"]])
model_msl3

model_msl6 <- summary(sim_aggregate6[["Model_MSL"]])
model_msl6

model_msl9 <- summary(sim_aggregate9[["Model_MSL"]])
model_msl9

model_msl12 <- summary(sim_aggregate12[["Model_MSL"]])
model_msl12

```

---

<!-- ## Optimal district magnitude -->

<!-- ```{r} -->
<!-- DM_min = ( log(0.05) - log(0.307) ) / -0.112 -->
<!-- DM_min -->

<!-- DM_max = ( log(0.1) - log(0.307) ) / -0.112 -->
<!-- DM_max -->

<!-- disp = 0.411 * exp(-0.080 * 9) -->
<!-- disp -->
<!-- ``` -->


# SEAT EXCESS AND SEAT BIAS

### 3-party system and Jefferson-d'Hondt

```{r}
# Data sampled using a truncated exponential distribution (the log-normal and uniform can also be used)
# Jefferson-d'Hondt
sb_dh3 <-
  simulate_Disp(
    seed = 0,
    dist = "lnorm", mean = 9, sd = 0.8,
    nd = 1,
    np = 3,
    ne = 500,
    formula = "dh",
    formula_dist = "hh",
    minTS = 3,
    maxTS = 20,
    jump = 1
  )

# Ranking
# subset(sb_dh3[["sb_bw"]], (elec %in% c("e1", "e11") ) & distTS == "12" )
```

----

```{r}
plot_dh3 <- plot_Disp(bias_data = sb_dh3, tse = c(0, 5/12, -1/12, -4/12))

plot_dh3[[1]]

# ggsave("sb_dh3_1.pdf", width = 22, height = 14, units = "cm")
```

----

```{r}
plot_dh3 <- plot_Disp(bias_data = sb_dh3, tse = c(0))

plot_dh3[[2]]

# ggsave("sb_dh3_2.pdf", width = 32, height = 12, units = "cm")
```

----

```{r}
plot_dh3[[3]]
```

----

```{r}
plot_dh3[[4]]
```

----

### 6-party system and Jefferson-d'Hondt

```{r}
# Data sampled using a truncated exponential distribution (the log-normal and uniform can also be used)
# Jefferson-d'Hondt
sb_dh6 <-
  simulate_Disp(
    seed = 0,
    dist = "lnorm", mean = 9, sd = 0.8,
    nd = 1,
    np = 6,
    ne = 10,
    formula = "dh",
    formula_dist = "hh",
    minTS = 3,
    maxTS = 20,
    jump = 1
  )

# Ranking
# subset(sb_dh6[["sb_bw"]], (elec %in% c("e1", "e2") ) & distTS == "3" )
```

----

```{r}
plot_dh6 <- plot_Disp(bias_data = sb_dh6, tse = c(0))

plot_dh6[[1]]

# ggsave("sb_dh6_1.pdf", width = 22, height = 14, units = "cm")
```

----

```{r}
plot_dh6 <- plot_Disp(bias_data = sb_dh6, tse = c(0))

plot_dh6[[2]]

# ggsave("sb_dh6_2.pdf", width = 32, height = 24, units = "cm")
```

----

```{r}
plot_dh6[[3]]
```

----

```{r}
plot_dh6[[4]]
```

----

### 9-party system and Jefferson-d'Hondt

```{r}
# Data sampled using a truncated exponential distribution (the log-normal and uniform can also be used)
# Jefferson-d'Hondt
sb_dh9 <-
  simulate_Disp(
    seed = 0,
    dist = "lnorm", mean = 9, sd = 0.8,
    nd = 1,
    np = 9,
    ne = 10,
    formula = "dh",
    formula_dist = "hh",
    minTS = 3,
    maxTS = 20,
    jump = 1
  )

# Ranking
subset(sb_dh9[["sb_bw"]], (elec %in% c("e1", "e2", "e10", "e11") ) & distTS == "3" )
```

----

```{r}
plot_dh9 <- plot_Disp(bias_data = sb_dh9, tse = c(0))

plot_dh9[[1]]

# ggsave("sb_dh9_1.pdf", width = 22, height = 14, units = "cm")
```

----

```{r}
plot_dh9 <- plot_Disp(bias_data = sb_dh9, tse = c(0))

plot_dh9[[2]]

# ggsave("sb_dh9_2.pdf", width = 32, height = 36, units = "cm")
```

----

```{r}
plot_dh9[[3]]
```

----

```{r}
plot_dh9[[4]]
```

----

### 3-party system and Sainte-Lague

```{r}
# Data sampled using a truncated exponential distribution (the log-normal and uniform can also be used)
# Sainte-Lague
sb_sl3 <-
  simulate_Disp(
    seed = 0,
    dist = "lnorm", mean = 9, sd = 0.8,
    nd = 1,
    np = 3,
    ne = 500,
    formula = "sl",
    formula_dist = "hh",
    minTS = 3,
    maxTS = 20,
    jump = 1
  )

# Ranking
# subset(sb_sl3[["sb_bw"]], (elec %in% c("e1", "e11") ) & distTS == "12" )
```

----

```{r}
plot_sl3 <- plot_Disp(bias_data = sb_sl3, tse = c(0))

plot_sl3[[1]]

# ggsave("sb_sl3_1.pdf", width = 22, height = 14, units = "cm")
```

----

```{r}
plot_sl3 <- plot_Disp(bias_data = sb_sl3, tse = c(0))

plot_sl3[[2]]

# ggsave("sb_sl3_2.pdf", width = 32, height = 12, units = "cm")
```

----

```{r}
plot_sl3[[3]]
```

----

```{r}
plot_sl3[[4]]
```

----

### 6-party system and Sainte-Lague

```{r}
# Data sampled using a truncated exponential distribution (the log-normal and uniform can also be used)
# Sainte-Lague
sb_sl6 <-
  simulate_Disp(
    seed = 0,
    dist = "lnorm", mean = 9, sd = 0.8,
    nd = 1,
    np = 6,
    ne = 10,
    formula = "sl",
    formula_dist = "hh",
    minTS = 3,
    maxTS = 20,
    jump = 1
  )

```

----

```{r}
plot_sl6 <- plot_Disp(bias_data = sb_sl6, tse = c(0))

plot_sl6[[1]]

# ggsave("sb_sl6_1.pdf", width = 22, height = 14, units = "cm")
```

----

```{r}
plot_sl6 <- plot_Disp(bias_data = sb_sl6, tse = c(0))

plot_sl6[[2]]

# ggsave("sb_sl6_2.pdf", width = 32, height = 24, units = "cm")
```

----

```{r}
plot_sl6[[3]]
```

----

```{r}
plot_sl6[[4]]
```

----

### 9-party system and Sainte-Lague

```{r}
# Data sampled using a truncated exponential distribution (the log-normal and uniform can also be used)
# Sainte-Lague
sb_sl9 <-
  simulate_Disp(
    seed = 0,
    dist = "lnorm", mean = 9, sd = 0.8,
    nd = 1,
    np = 9,
    ne = 10,
    formula = "sl",
    formula_dist = "hh",
    minTS = 3,
    maxTS = 20,
    jump = 1
  )

```

----

```{r}
plot_sl9 <- plot_Disp(bias_data = sb_sl9, tse = c(0))

plot_sl9[[1]]

# ggsave("sb_sl9_1.pdf", width = 22, height = 14, units = "cm")
```

----

```{r}
plot_sl9 <- plot_Disp(bias_data = sb_sl9, tse = c(0))

plot_sl9[[2]]

# ggsave("sb_sl9_2.pdf", width = 32, height = 36, units = "cm")
```

----

```{r}
plot_sl9[[3]]
```

----

```{r}
plot_sl9[[4]]
```

----

### 3-party system and Imperiali

```{r}
# Data sampled using a truncated exponential distribution (the log-normal and uniform can also be used)
# Imperiali
sb_imp3 <-
  simulate_Disp(
    seed = 0,
    dist = "lnorm", mean = 9, sd = 0.8,
    nd = 1,
    np = 3,
    ne = 10,
    formula = "imperiali",
    formula_dist = "hh",
    minTS = 3,
    maxTS = 20,
    jump = 1
  )

```

----

```{r}
plot_imp3 <- plot_Disp(bias_data = sb_imp3, tse = c(0))

plot_imp3[[1]]

# ggsave("sb_imp3_1.pdf", width = 22, height = 14, units = "cm")
```

----

```{r}
plot_imp3[[2]]

# ggsave("sb_imp3_2.pdf", width = 32, height = 36, units = "cm")
```

----

```{r}
plot_imp3[[3]]
```

----

```{r}
plot_imp3[[4]]
```

----

### 6-party system and Imperiali

```{r}
# Data sampled using a truncated exponential distribution (the log-normal and uniform can also be used)
# Imperiali
sb_imp6 <-
  simulate_Disp(
    seed = 0,
    dist = "lnorm", mean = 9, sd = 0.8,
    nd = 1,
    np = 6,
    ne = 10,
    formula = "imperiali",
    formula_dist = "hh",
    minTS = 3,
    maxTS = 20,
    jump = 1
  )

```

----

```{r}
plot_imp6 <- plot_Disp(bias_data = sb_imp6, tse = c(0))

plot_imp6[[1]]

# ggsave("sb_imp6_1.pdf", width = 22, height = 14, units = "cm")
```

----

```{r}
plot_imp6[[2]]

# ggsave("sb_imp6_2.pdf", width = 32, height = 36, units = "cm")
```

----

```{r}
plot_imp6[[3]]
```

----

```{r}
plot_imp6[[4]]
```

----

### 9-party system and Imperiali

```{r}
# Data sampled using a truncated exponential distribution (the log-normal and uniform can also be used)
# Imperiali
sb_imp9 <-
  simulate_Disp(
    seed = 0,
    dist = "lnorm", mean = 9, sd = 0.8,
    nd = 1,
    np = 9,
    ne = 10,
    formula = "imperiali",
    formula_dist = "hh",
    minTS = 3,
    maxTS = 20,
    jump = 1
  )

```

----

```{r}
plot_imp9 <- plot_Disp(bias_data = sb_imp9, tse = c(0))

plot_imp9[[1]]

# ggsave("sb_imp9_1.pdf", width = 22, height = 14, units = "cm")
```

----

```{r}
plot_imp9[[2]]

# ggsave("sb_imp9_2.pdf", width = 32, height = 36, units = "cm")
```

----

```{r}
plot_imp9[[3]]
```

----

```{r}
plot_imp9[[4]]
```
